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1 . INTRODUCTION 


The use of ceramics in heavy duty diesel engine applications is 
especially promising [1]. It has been shown in the TACOM/Cummins Adiabatic 
Engine Program that reductions in heat losses at appropriate points in the 
diesel engine system result in substantially increased exhaust enthalpy [2]. 
One of the most promising concepts for taking advantage of this increased 
exhaust enthalpy is the turbocharged turbo-compounded diesel engine cycle [3]. 

This engine concept consists of several sub-systems: compressor, engine, 

turbocharger turbine, compounded turbine, ducting and heat exchangers. 
Applications of ceramic materials to different sub-system components have 
widely different degrees of difficulty. Therefore, there is a need for a 
computer simulation of this engine concept, at the appropriate level of 
detail, to enable engineers to define the trade-offs associated with 
introducing ceramic materials in various parts of the total engine system, and 
to carry out system optimization studies. 

This report describes the model being developed for simulating the 
behavior of the total engine system. The focus of this total system 
simulation is to define the mass and energy transfers (heat transfers, heat 
release, work transfers) in each sub-system and the relationship between the 
sub-systems. Since this system model must function as a single-unit, in this 
development program, a deliberate effort is made to maintain a balance in the 
complexity of the various sub-system descriptions. 

Figure 1 illustrates the proposed overall model structure. The air flow 
is followed through an air filter, ducting, turbocharger compressor, ducting, 
cooler and engine intake system, to the diesel engine. The engine operating 
cycle is then modelled. The exhaust gas flow is followed through the engine 
exhaust ports, manifold and ducting to the turbocharger turbine; then followed 
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through additional ducting to the compounded turbine, the exhaust system and 
muffler to the atmosphere. Engine friction sub-models are then used to obtain 
brake quantities from the computed indicated quantities. 

The engine model of the total system simulation is a mathematical model 
of the processes occur ing in a diesel engine over a complete cycle. The 
direct-injection four-stroke diesel-engine cycle-simulation developed in this 
work is based on two existing M.I.T. engine cycle simulations: (i) a spark- 

ignition engine simulation [4], and (ii) a divided-chamber diesel engine 
simulation [5]. The direct-injection diesel simulation has retained many of 
the basic modelling assumptions (gas exchange process, thermodynamic and 
transport properties routines), as well as the modular structure and efficient 
computational algorithm of the parent codes. 

The tasks in the first year of the program were to develop this DI diesel 
simulation, modify this single cylinder engine model to predict the 
performance of a multi-cylinder engine, and to implement the appropriate logic 
for coupling the multi-cylinder engine with the turbocharger, compound turbine 
and other system component models. 

The purpose of this report is to summarize the basic assumptions and 
mathematical relationships used in this simulation of the turbocharged 
turbo-compounded diesel engine cycle. The report is arranged as follows. 
Sections 2, 4 and 5 describe the simulation of the diesel engine itself. 
Sections 6 and 7 describe the models used for the other system components: 
the turbocharger compressor, and turbine, compound turbine, wastegate, cooler, 
manifold and ducting nodels. Sections 3, 8 and 9 relate to all parts of the 
total engine system. Sections 10 and 11 describe the integration together of 
all the subsystems, and program inputs and outputs. 
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2. BASIC ASSUMPTIONS OF ENGINE MODEL 


The diesel engine cycle is treated as a sequence of continuous processes: 
intake, compression, combustion (including expansion), and exhaust. The 
duration of the individual processes are as follows. The intake process 
begins when the intake valve opens (IVO) and ends when the intake valve closes 
(IVC). The compression process begins at IVC and ends at the time of ignition 
(IGN). The combustion process begins when ignition occurs and ends when the 
exhaust valve opens (EVO). The exhaust process begins at EVO and ends at IVO 
(and not when the exhaust valve closes). 

In the engine simulation, the system of interest consists of the 
instantaneous contents of a cylinder, i.e. air and combustion products. In 
general, this system is open to the transfer of mass, enthalpy, and energy in 
the form of work and heat. Throughout the cycle, the cylinder is treated as a 
variable volume plenum, spatially uniform in pressure. Furthermore, the 
cylinder contents are represented as one continuous medium by defining an 
average equivalence ratio and temperature in the cylinder at all times. 

Gas properties are obtained assuming ideal gas behavior. At low 
temperatures (below 1000 K), the cylinder contents are treated as a 
homogeneous mixture of non-reacting ideal gases. At high temperatures (above 
1000 K), the properties of the cylinder contents are calculated with allowance 
for chemical dissociation by assuming that the burned gases are in 
equilibrium, using an approximate calculation method based on hydrocarbon air 
combustion. 

Quasi-steady, adiabatic, one-dimensional flow equations are used to 
predict mass flows past the valves. The intake and exhaust manifolds are 
treated as plenums whose pressure and temperature histories are determined by 
solution of the manifold state equations. When reverse flow past the intake 
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valve occurs, a rapid mixing model tor the back flow gases is assumed. The 
engine model couples with the compressor and heat exchanger at the inlet 
manifold, and with the turbocharger turbine at the exhaust manifold. 

The compression process is defined so as to include the ignition delay 
period; i.e., the time interval between the start of the injection process 
(the point at which the injector needle starts to lift) and the ignition 
point (the start of positive heat release due to combustion). The total 
length of the ignition delay is related to the mean cylinder gas temperature 
and pressure during the delay period by an empirical Arrhenius expression. 

The combustion process is modelled as a uniformly distributed heat 
release process. The rate of heat release is assumed to be proportional to 
the rate of fuel burning which is modelled empirically. Since the diesel 
combustion process is comprised of a pre-mixed and diffusion-controlled 
combustion mechanism, it was decided to use Watsop's fuel burning rate 
correlation [6], consisting of the sum of two algebraic functions, one for 
each combustion mechanism. The proportion of the total fuel injected that 
is burnt by either mechanism depends on the length of the ignition delay 
period and the engine load and speed. 

Heat transfer is included in all the engine processes. Convective heat 
transfer is modelled using available correlations for turbulent flow in 
pipes. The characteristic velocity and length scales required to evaluate 
these correlations are obtained from a mean and turbulent kinetic energy 
model. Radiative heat transfer is added during combustion. 

Finally, an empirical friction model is used to convert the indicated 
engine performance quantities to brake performance quantities. 
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3. CONSERVATION EQUATIONS 


In this section, equations for the conservation of mass and energy will 
be stated for the contents of a generalized control volume treated as an open 
system. The conservation equation for the fuel mass will be used to develop a 
differential equation for the change in fuel-air equivalence ratio of the 
system. The energy conservation equation will be expanded to obtain a 
differential equation for the change in temperature of the thermodynamic 
system. 

3.1 Conservation of Mass 

3.1.1 Conservation of total mass 

The rate of change of the total mass in any control volume or open system 
is equal to the sum of the mass flow rates into and out of the control volume, 

i.e. 

m = I m (3-1) 

Note that the convention used in our model assumes that mass flow rates 
into the control volume are taken as positive, while mass flow rates out of 
the control volume are negative. 

3.1.2 Conservation of fuel mass 

In particular, conservation of the fuel species can be expressed as 

• • 

m = Z m (3-2) 

j 

where m denotes the fuel content in the control volume (includes fuel added 
by injection and fuel in the form of combustion products). 

Defining the fuel fraction, F, in the control volume as 

F = m f /m 
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where m is the total mass in the control volume, equation (3~2) can be 
re-written as 

3j(*F) - t i jFj (3-3) 

where Fj denotes the fuel fraction of the mass flow entering or leaving the 
control volume through the j port. 

Differentiating the left-hand side of equation (3-3) gives 

mF = E m.F. - mF ( 3 -H) 

. J J 
j J J 

or substituting for m from equation (3~ 1 ) results in a differential equation 
for the change in the fuel fraction of the control volume, i.e. 



in in 

i 3i 

^ = FASTO 

where m is the mass of air in the control volume and FASTO denotes the 
a 

stoichiometric fuel to air ratio. 

Expressing the equivalence ratio, $, in terms of the fuel fraction, F, 

i.e. 


1 F 

* = FASTO T-F (3_7) 

and differentiating ( 3~ 7 ) with respect to time we obtain an equation for the 
rate of change of the equivalence ratio of the control volume, i.e. 
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' 1 F 

* = 2 (3-8) 

FASTO (1-F) 

with F given from equation (3'5) . 

3.2 Conservation of Energy 

The general energy equation for an open thermodynamic system may be 
written as 


E = E m.h Q - W (3-9) 

j J J 

with the rate of change of the energy of the system being given by 

E - SjjCnt.) - |(pV) (3-10) 

where 

E m.h. is the net rate of influx of enthalpy 
j J J 

Q = E Q. is the total heat transfer to the walls, i.e. the sum of the 

W . 1 

1 

heat transfer rates to the different surfaces of the control 
volume of interest 

W = pV is the rate at which the system does work by boundary 
displacement. 

The dots denote differentiation with respect to time. Note that the 
convention used is that heat loss from the system and work done by the system 
are taken as positive. 

Differentiating the left-hand side of equation (3~9) gives 
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mh = E ffl.h - Q + pV - mh (3-11) 

J J j -w 

The contents of any control volume, i.e. air and combustion products, can 
be represented as one continuous medium by defining an average equivalence 
ratio at all times. Gas properties are obtained assuming ideal gas behavior 
and thermodynamic equilibrium. With these assumptions, we can express the 
enthalpy, h, and the density, p, of the mixture of air and combustion products 
as 

h = h(T,p, <*>) (3-12) 

p = p(T,p,<|>) (3-13) 

Hence, the rate of change of the above fluid properties with respect to time, 
or crank-angle, can be written as 


where 


n - c p t . c T p * c ( 


, 9h , 

c = ( 7 - ) , 

P 3T p , <j 

V dp i, <J) 

,9h. 

V T,P 


and 


P = A T + (InV A P + ^SaV n $ 

3T p , 9p T , 4> 94> T,p 

The equation of state for ideal gases 
p = RpT 

can be expressed in differential form as 


(3-14) 


(3-15) 


(3-16) 
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2 = 5 + £ + 1 

p R p T 

Re-arranging equation (3-17) and using (3~16), we can write 


(3-17) 


R = nT P ~ ~p~ p " -B p T 
PT P T PT 2 


Substituting for p from equation ( 3“ 1 5 ) into the above equation, 


(3-18) 


we can 


express R in terms of p, T and <J>, i.e. 


r = (1 2- 9fi)n - (-2- h- -2- 3% - -2_ iG i 

1,1 pi 2 p 2 ! 31 p 2 T 3 * 


(3-19) 


From the differential form of tine equation of state, we can express the 
time rate of change of pressure as 


,R m T V, 
p P R + m + T V ) 


( 3 - 20 ) 


or substituting for R from equation (3~19), and with some manipulation, we 
obtain the time rate of change of pressure: 


n = P r- - I 22 r - I 1 + Si 

p - 3p/3p V p 9T T p 34- * m 


(3-21 ) 


Returning to the energy equation (3 - 11), and expressing h in terms of 
its partial derivatives with respect to T, p and <p, we obtain 


mc_T = E m.h.- Q + (V - mc„)p - me, <p - hm 

Pjjjw 


( 3 - 22 ) 


Finally, substituting for p from equation (3-21), we obtain an 
equation for the rate of change of temperature that does not explicitly depend 
on the rate of change of pressure of the system, i.e. 
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me T = (E m.h,- Q ) - me <|> - hm + 
P j J J w 4> y 


* -X- If - 1 f* T - 1 HI . S, 
V p dT p d<(> m 


“ here B ■ rs'p/spj 0 - “V 


or dividing (3-23) by m, and collecting terms in T, $ and m get 


(c + ; ||)T = - (E m.h.- Q ) - (c + - | fi ) j, + 
P P 3T m . jj \ <J> p 3f T 


+ B - (1 
m 


h) _ o Y 

B 8 V 


leading to 


* = K (1 - §> - hi ♦ + k (z ; j h j- w 


where 


+ § ifi , c + LMm.fl . c 

P 3T p Op/ap) v p C T 

C + 8 ^ = c + LMML1 _ c 

<{> P 3<j> 4> (3p/3p) P T 



(3-23) 

(3-24) 


(3-25) 


(3-26) 


(3-27) 

(3-28) 


10 



i&INAL PAGE IS 
: POOR QUALITY 


. APPLICATION OF CONSERVATION EQUATIONS TO ENGINE PROCESSES 


The conservation equations derived in section 3 can be applied to the 
four engine processes as follows: 
i) Intake 


m = in. - m 
m ex 


• m . m 

F = (F. - F) - ---(F - F) 

ram ra ex 


T -f[5(1 -§) -l-H + U*. h. -m h -0 )] 

Am B V B Bra ■ in m ex ex w 


ii) Compression 


m = 0 


F = 0 


* R V Q w 

T = -[ -] 

A V Bm J 


(iii) Combustion 


F = ~(1 - F) 
ra 


* - aV (1 - s> - hi ♦ + k ( ^f- V ] 


where h^, is the absolute fuel enthalpy. 


n 




(iv) Exhaust 


port 

flow 



F = - ---(F - F) 
m ex 


1 B r ex,, h> V C ' 1 , * : 

T = -[- (1 “5) _ v _ n < J ) + m h " Q )] 

A m B V B Bm ex ex w 


Note that the fuel fraction of the mass flow rate through the 
F ex> could be different from the cylinder fuel fraction, F, 
situations. 


( 4 - 10 ) 

( 4 - 11 ) 

( 4 - 12 ) 

exhaust 
in reverse 
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5 . MODELLING OF ENGINE PROCESSES 


5.1 Gas Exchange 

A one-d imens ional quasi-steady compressible flow model is used to 
calculate the mass flow rates through the intake and the exhaust valves during 
the gas exchange process. The manifolds are treated as plenums with known 
pressures. Furthermore, the temperature and average equivalence ratio of the 

l 

intake charge (fresh air and EGR at intake manifold conditions) and the 
exhaust charge (mixture of air and combustion products at cylinder conditions) 
are known. When reverse flow into the intake manifold occurs, a rapid-mixing 
model is used, i.e. perfect and instantaneous mixing between the back-flowing 
charge and the intake charge is assumed. 

At each step of the gas exchange process, values for the valve open areas 
and discharge coefficients are obtained from tabulated data. Given the open 
area, the discharge coefficient, and the pressure ratio across a particular 
valve, the mass flow rate across that valve is calculated from: 


P„ , , e, 2/T P, <Y*1)/Y 1/2 

" - V 5F^?V-Tt<- p S > - if) 3) 

o o o 


where 


(5-D 


c^ = discharge coefficient 
A = valve open area 

P q = stagnation pressure upstream of restriction 
P s = static pressure at restriction 
T q = stagnation temperature upstream of restriction 
Y = ratio of specific heats 
R = gas constant 

When the kinetic energy in the cylinder is negligible, the stagnation 
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pressure and temperature reduce to the static pressure and temperature, 
respectively. 

For the case of choked flow, equation (5-1) reduces to 
p (Y+1 )/2(Y-1) 

m - ) (5-2) 

O 

The mass, m(t), in the cylinder at any time t can he found from 
integration of the mass conservation equation (3-1), i.e. 

f * t * 

m(t) = m + / m. dt - / m dt (5-3) 

o t in t ex 
o o 

where m Q is the mass in the cylinder at cycle start (IVO). 

5.2 Combustion Model 

The diesel combustion process is a very complex, unsteady, heterogeneous, 
three-dimensional process. A good mathematical combustion analysis would 
require accurate models of compressible viscous air motion, fuel spray 
penetration, droplet break-up and evaporation, air entrainment, combustion 
kinetics, turbulent diffusion and so on. The details of the combustion 
process would depend on the characteristics of the fuel, the design of the 
combustion chamber and the fuel injection system, and on the engine's 
operating conditions. Although an adequate conceptual understanding of diesel 
combustion has been developed to date, a comprehensive quantitative model of 
all the individual processes has not yet been proposed. 

A relatively successful approach to the problem of combustion simulation 
has been to model combustion as a heat release process, as originally proposed 
by Lyn [7]. The rate of heat release (or, equivalently, the rate of fuel 
burning) can be defined as the rate at which the chemical energy of the fuel 
is released by combustion. Based on heat release analysis and high-speed 
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photography, Lyn has provided an excellent description of the different stages 
of diesel combustion. These stages can be identified on the typical rate of 
heat release diagram for a DI engine shown in Fig. 2 as follows: 

Ignition delay (ab): The period between the dynamic injection point (the 

point at which the injector needle starts to lift) and the ignition point (the 
start of positive heat release due to combustion). 

Pre-mixed or rapid combustion phase (be): In this phase, combustion of the 

fuel which has mixed with air to within the flamability limits during the 
ignition delay period occurs rapidly in a few crank angle degrees. This 
results in the high initial rate of burning generally observed in direct- 
injection diesel engines. 

Mixing controlled combustion phase (cd): Once the fuel and air premixed 

during the ignition delay have been consumed, the heat release rate (or 
burning rate) is controlled by the rate at which mixture becomes available for 
burning. The heat release rate may or may not reach a second (usually lower) 
peak in this phase; it decreases as this phase progresses. 

Late combustion phase (de): Heat release continues at a low rate well into 

the expansion stroke. Eventually, the burning rate asymptotically approaches 
zero. The nature of combustion during this phase is not well understood. 
Possible processes are that a snail fraction of the fuel may not yet have 
burned, or energy present in soot and fuel-rich combustion products can still 
be released, etc. Given the somewhat arbitrary limits of this phase, 
combustion models usually focus on the main heat release periods, i.e. the 
pre-mixed and mixing controlled combustion phases. 
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Based on his observations, Lyn [7] proposed a method of predicting 
burning rates from the rate of fuel injection. The fuel injected is divided 
into elements according to the order in which they enter the combustion 
chamber. These fuel elements become "ready for burning" according to a 
certain law (with increasing burning time as injection proceeds). Thus, a 
"ready for burning" diagram can be obtained from the rate of injection 
diagram. At the ignition point, which occurs after the lapse of the delay 
period, the part of the injected fuel which has been made "ready for burning" 
is added on the total "ready for burning" diagram, causing the sharp peak in 
the burning rate diagram. Subsequent burning is essentially proportional to 
the rate of injection. Although this method gives a reasonable fit to data 
obtained over a wide range of speeds and loads, it requires further refinement 
and calibration before it can be used in computer simulation work. 

It is necessary in combustion modelling, in the context of computer 
simulation to predict the performance of new engine concepts, to describe the 
apparent fuel burning rate by algebraic expressions. The constants in these 
expressions can be chosen suitably to reflect the dependence of the actual 
fuel burning rate on engine type and particular operating conditions. 

Ship inski et al [8] attempted to correlate the apparent rate of fuel 
burning with the rate of fuel injection by fitting a Wiebe function to heat 
release diagrams obtained from tests on a high-speed swirl-type direct 
injection engine. Although Shipinski obtained a reasonable agreement with his 
engine data, the heat release shape defined by the Wiebe function alone has a 
notable difference from the two-part characteristic with the initial spike 
that is measured on most engine types. 

To overcome this problem, Watson et al [9] proposed that the apparent 
fuel burning rate could be expressed as the sum of two components, one 
relating to pre-mixed and the other to diffusion-controlled burning, i.e. 
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(5-4) 


where m is the fuel burning rate with respect to crank angle 

and subscripts t, p, d denote total, pre-mixed and diffusion burning, 

respectively. 

In order to quantify the proportion of the fuel burnt by either 
mechanism, a phase proportionality factor, 6, is introduced. This expresses 
the cumulative fuel burnt by pre-mixed burning as a fraction of the total fuel 
injected, i.e. 


m 

B = ~ E (5-5) 

m t 

Consequently, a non-dimensional apparent fuel burning rate curve can be 

\ 

written as 

M fc (t) = BM (t) + (1 - B)M (x) (5-6) 

o p u 

where M(x) = non-dimensional burning rate distribution 
and x = normalized crank position. 

The phase proportionality factor, S, is considered to be controlled by 
the length of the ignition delay period (since the fuel that is prepared for 
burning during this period governs the pre-mixed burning phase), and the 
overall cylinder equivalence ratio, <j> . This can be expressed by a 

relation of the form 


4 > 

8=1- a---- 

ID° 


where ID = ignition delay (see Section 5.3)* 


(5-7) 
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From correlation with data from turbo-charged truck engines, Watson [9] 
obtained the following range of values for a, b, and c: 

0.8 < a < 0.95 

0.25 < b < 0.45 (5-8) 

0.25 < c < 0.50 

Futhermore, Watson concluded that the best representation of the 
experimental data was achieved using the following component burning rate 
distributions: 

Premixed burning: 


M (t) = 1 - (1 - t p1 ) Cp2 


(CrD c , (c_-i ) 


or M ( t) = C .C _t 
P Pi p2 


Pi 


(1 - , P1 ) P2 


(5-9a) 
(5— 9b ) 


Diffusion controlled burning (Wiebe function) 

C 

M d (i) = 1 - exp(-C d1 T ) (5-1 0a) 

(c -i) C 

or ( t ) = C d1 C d2 T exp(-C d1 T ) (5-1 0b) 

where C ,, C 0 , C , C.„ are shape factors, 
pi p2 di d2 

The shape factors in equations (5-9) and (5-10) can be determined as a 
function of the engine operating conditions. Using data from a typical truck 
engine, Watson established the following correlations: 


C pl = 2.0 + 1 .25 x 10 8 (ID x N) 2,4 

(5-1 la) 

C p2= 5000 

(5-1 1b) 

C = 14.2 4> "O- 6 ^ 4 

di ove 

(5-1 1c) 
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0.25 


(5-1 Id) 


C d2= 


0.79 C 


dl 


where 

ID = ignition delay (see Section 5.3) 

N = engine speed (RPM) 

<j>ov e = overall cylinder equivalence ratio. 


Watson's correlation gives good agreement between experimental and 
predicted results over the entire range of engine operating conditions, and is 
flexible enough to be applied to various types of diesel engines. Thus, we 
have followed Watson's approach to describe the heat release profile for our 
engine simulation. If necessary, the constants (or the form) of correlations 
(5-7) and (5-11) will be modified to fit experimental data for our engine 
design. 


5.3 Ignition Delay Model 

The ignition delay period in a diesel engine was defined in section 5.2 
as the time (or crank angle) interval between the start of injection and the 
start of combustion. The. start of injection is usually taken as the time when 
the injector needle lifts off its seat (determined from a needle lift 
indicator). The start of combustion is more difficult to determine. It is 
best identified from the change in slope of the heat release rate which occurs 
at ignition. 

Both physical and chemical processes must take place before the injected 
fuel can burn. The physical processes are: the atomization of the liquid 

fuel jet; the vaporization of the fuel droplets; the mixing of fuel vapor with 
air. The chemical processes are the preccrabust ion reactions of the fuel, air, 

{ 
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residual gas mixture which lead to autoignition. These processes are affected 
by engine design and operating variables, and fuel characteristics. 

Ignition delay data from fundamental experiments in combustion bombs and 
flow reactors have usually been correlated by equations of the form 

ID = A p~ n exp(E/RT) (5-12) 

where ID = ignition delay [ms] 

E = apparent activation energy for the fuel autoignition process 
R = universal gas constant 
p = gas pressure [atm] 

T = gas temperature [K] 

and A and n are constants dependent on the fuel and the injection and 
air-flow characteristics. Representative values for A, n and E are given in 
Table 1 . 

These correlations have usually been derived from tests in uniform air 
environments where the pressure and temperature only changed due to the 
cooling effect of the fuel vaporization and fuel heating processes. However, 
in a diesel engine, pressure and temperature change considerably during the 
delay period due to the compression resulting from piston motion. .Another 
problem is that the form of these simple correlations is not sufficiently 
flexible to allow all the influencing fuel and engine parameters to be 
included in the calculation of the ignition delay. 

Hardenberg and Hase [10] have developed an empirical formula for 
predicting the ignition delay in DI engines as a function of fuel 
characteristics, engine parameters and ambient conditions. Dent [11] has 
shown that this formula can give reasonable agreement with experimental data 
over a wide range of engine conditions. However, the pressure and temperature 
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used in this correlation are identified as the corresponding conditions at top 
deal centre, estimated using a polytropic model for the compression process. 

It is felt that such a polytropic model is not appropriate in our simulation 
context, where pressures and temperatures can be accurately predicted over the 
duration of the ignition delay period. Therefore, we propose to use a formula 
similar to Equation (5-12), with pressure and temperature taken at their 
arithmetic mean values during the delay period to account for changing 
conditions. 


5.4 Heat Transfer 

The heat transfer mechanisms in a diesel engine include forced convection 

(Q ) from the turbulent flow in the cylinder to the combustion chamber walls, 

and radiation (Q ) from the flame and the burning soot particles. The total 

heat transfer rate (Q ) is therefore given by 

w 




(5-13) 


The convective heat transfer at the gas-to-cylinder wall interface will 
depend on the temperature gradient in the boundary layer at the surface. 
However, due to the inherent difficulties in calculating the boundary layers 
in a diesel engine, the convective heat transfer rate is usually expressed as 


Q = h A (T - T ) (5-14) 

c gw 

where 

h = : convective heat transfer coefficient 
A = surface area 
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bulk mean gas temperature 


T = 
g 

T w = local surface temperature on cylinder wall, head, or piston, as 
appropriate. 

The problem is then to devise a method to calculate the convective heat 
transfer coefficient that appears in (5-1 *1). The approach usually taken is to 
calculate h from a Nusselt-Reynolds number correlation analagous to that used 
for steady turbulent flow in a pipe [12-1 6], i.e. 

d e 

Nu = a Re Pr (5-15) 

where 

Nu = hL/A : Nusselt number 

Re = VL/v : Reynolds number 

Pr = pc /A: Prandtl number 
P 

L = a characteristic length 
V = a characteristic velocity 
A = thermal conductivity 
v = kinematic viscosity 
y = dynamic viscosity 
p = density 

0 ^ = specific heat at constant pressure 
and a, d, e are constants adjusted to fit experimental data. 

Fortunately, there is little variation in the Prandtl number for air and 
combustion products, which is usually close to unity. Consequently, we may 
drop the Prandtl number dependence in equation (5-15) with little loss in 
accuracy, so that 

Nu = a Re d (5-16) 

To calculate the convective heat transfer coefficient from correlation 
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(5-16), instantaneous values for the characteristic length and velocity 
scales, and the gas transport properties ( pi, p, and X) are needed. Currently, 
it is not possible to predict these parameters and their spatial and temporal 
variation with any accuracy in an internal combustion engine. To overcome 
this difficulty, representative values of the characteristic length and 
velocity scales, and the gas temperature, pressure, and equivalence ratio at 
which the gas properties are to be evaluated are chosen. 

For our heat transfer model, the characteristic length scale is taken to 
be the macroscale of turbulence, as defined by equation (5-28). (See 
turbulent flow model) . The characteristic velocity V is postulated to be an 
effective velocity due to contributions from the mean kinetic energy, the 
turbulent kinetic energy and piston motion, i.e. 



a = 0.035 to 0.13 (5 _ 18 

d = 0.7 to 0.8 

depending on intensity of change motion. The gas density and the transport 
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properties, u and A, that appear in correlation (5-16) are evaluated at the 
mean gas temperature, pressure, and equivalence ratio. 

Radiative heat flux is significant in diesel engines. However, due to 
the difficulty of measuring instantaneous flame temperature and heat flux, it 
is not easy to quantify the contribution of radiation to the overall heat 
transferred. Estimates of the relative importance of radiation have varied 
between a few and 30 percent of total heat transfer, and vary according to 
engine type. [13,14,17,18]. 

Annand [13] has expressed the radiation term as 


4 4 

Q = k (T - T ) 
r r g w 


(5-19) 


where 

= empirical radiation constant 

T = bulk mean gas temperature 
8 

T w = local surface temperature on cylinder wall, head or piston, 
as appropriate 

During the intake, compression and exhaust processes, when the radiative 
heat flux would be zero, k^ = 0. During combustion, .Annand and Ma [14] 
suggested that 

k = C o 
r r 

-12 2 4 

where o = Stephan-Boltzmann constant (56.7x10 kw/m K ), 
and is an adjustable constant with values in the range of 1.3 - 3.1, 

depending on the engine speed and load. The value of C^, exceeding unity, 
suggests that the actual radiation temperature is appreciably above the bulk 
mean value. 

Spatial variations in heat transfer rates are accounted for by assigning 
different surface temperatures to the cylinder walls, head, piston, etc., 
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based on practical experience. If these temperatures are to be calculated, a 
heat conduction model is required for the piston, valves, cylinder head, and 
liner. 

The total instantaneous heat tranfer rate will be calculated from the sum 
of the convective terms (predicted by equation 5-14) plus the sum of the 
radiative terms (predicted by equation 5-19). 

5.5 Turbulent Flow Model 

The heat transfer model of the cycle simulation requires estimates of the 
characteristic velocity and length scales. To estimate these scales in a way 
which incorporates the key physical mechanisms affecting charge motion in the 
cylinder, a turbulent flow model is used. This model is a variation of the 
models used by Manscuri [5] and Poulos [4] in previous engine simulation work. 

The turbulence model consists of a zero-dimensional energy cascade. Mean 
flow kinetic energy, K, is supplied to the cylinder through the valves. Mean 
kinetic energy, K is converted to turbulent kinetic energy, k, through a 
turbulent dissipation process. Turbulent kinetic energy is converted to heat 
through viscous dissipation. When mass flows out of the cylinder, it carries 
with it both mean and turbulent kinetic energy. Figure 3 illustrates the 
energy cascade model. 

At any tin® during the cycle, the mean flow velocity, U, and the 
turbulent intensity, u', are found from knowledge of the mean and turbulent 
kinetic energies, K and k, respectively. Thus, the following equations apply: 

K = | m U 2 (5-20) 

k = | m u’ 2 (5-21) 
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where the factor 3 in equation (5-21) comes from assuming that the small scale 
turbulence is isotropic (and accounting for all three orthogonal fluctuating 
velocity components). 

The time rate of change of the mean kinetic energy, K is given by 


lm.V. 2 - P -K-§ 

dt 2 1 1 m 


( 5 - 22 ) 


Similarly, the rate of change of the turbulent kinetic energy, k, is 


with 

where 


— = P - me - k -- + A 
dt m 


e = 


u' 3 (2k/3m) 372 


l l 

m = mass in the cylinder 


m i = mass flow rate into the cylinder 


(5-23) 

(5-24) 


mass flow rate out of the cylinder 

V^= jet velocity into the cylinder 

P = rate of turbulent kinetic energy production 

e = rate of turbulent kinetic energy dissipation per unit mass 

A = rate of turbulent kinetic energy amplification due to rapid 

distortion 

l = characteristic size of large-scale eddies 
In equations (5-22) and (5-^3). the production term P has to be defined 
in terms of flow and geometrical parameters of the chamber. However, since 
the above model does not predict spatially resolved flow parameters, P must be 
estimated from mean flow quantities only. 
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Assuming that turbulence production in the engine cylinder is similar to 
turbulence production in a boundary layer over a flat plate [19], we can 
express P as 


- ■ 


where u = C k /(me) is turbulent viscosity, 
« u 


and C = 0.09 is a universal constant, 

u 


(5-25) 


Again, as the velocity field in the cylinder is not known, the velocity 
gradient (3U/3y) is approximated as 


<38, 2 . c <-5 2 


(5-26) 


•nere C is an adjustable constant and L is a geometric length scale. 

P 

Using equations (5-25), (5-26), (5-20) and (5-24), we can express P as 


3/2 1 /2 
7 ■ X? k 

p = 2(§) c c fl ra(£r 

2 u g ^2 m 


(5-27) 


Furthermore, the characteristic size of the large-scale eddies, l, and 
the representative geometric length scale, Lwill be both identified with the 
macroscale of turbulence, assumed to be given by 

i = L = V/(ttB 2 /4) (5-28a) 

where V is the instantaneous volume of the combustion chamber and 3 is the 
cylinder bore, subject to the restriction that 

L < B/2 f 

Hence, equation (5 _ 27) can be re-written as 
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(5-28b) 



During the compression and the combustion processes, the turbulent 
kinetic energy decays due to viscous dissipation. On the other hand, at the 
same time the turbulent kinetic energy is amplified due to the rapid 
distortion that the cylinder charge undergoes with rising cylinder pressures. 
Consequently, an amplification term, A, was added to equation (5-23) to 
account for this effect. The amplification term will be larger during 
combustion when the unburned gas is assumed to be compressed by the flame at a 
sufficiently high rate. However, in the diesel engine context (high 
compression ratios), it is proposed to keep the amplification term during 
compresson, too.. 

Using equation (5-21), the rate of turbulence amplification due to rapid 
distortion can be expressed as 


A = 3 m u 



(5-30) 


where the rate of change of the turbulent intensity du'/dt can be estimated 

assuming that conservation of mass and angular momentum can be applied to the 

large scale eddies during the rapid distortion period. 

Under these assumptions, conservation of mass for a single eddy of volume 

V requires that 
L 


p V = p V T 
L o Lo 


(5-3D 


where p is the mean gas density and subscript o refers to the conditions at 
the start of compression. Then, since 


V 


L 


(5-32) 


where L is the macroscale of turbulence, we can re-write (5-31) as 



(5-33) 


1/3 


~ = (--) 
L o p 


Conservation of eddy angular momentum requires that 

U L = U L 
(l) 0)0 o 

where is the characteristic velocity due to eddy vorticity. 

Combining (5-33) and (5-34) with the assumption that 

U - u* 

0) 

the following relation is obtained for the evolution of the turbulent 
intensity during the rapid compression period 

, 1/3 

, (fL) 

u' p 
o o 

Differentiating both sides of equation (5-36) and re-arranging we get 


(5-34) 


(5-35) 


(5-36) 


du'_ u]_ dg 
dt 3p dt 


(5-37) 

Hence, combining equations (5~30) and (5-37), the rate of turbulence 


amplification is given by 


. _ mu^_ dg 
' p dt 

or, introducing (5-21), 


(5-38) 


A = | k 
3 p 


(5-39) 


with p given from equation (3~15). 
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5.6 Engine Friction Model 


To convert indicated engine performance quantities to brake engine 
performance quantities, engine friction calculations are required. However, 
the measurement and analysis of engine frictional losses has never been 
satisfactorily resolved. This is primarily due to the inherent problem of 
direct, accurate measurement of these losses under actual running conditions. 
This problem occurs because the total loss is a summation of losses arising 
from the operation of the many components of the engine, and these components 
respond differently to changes in pressure, temperature and speed. 

Direct motoring of an engine is the common method of measuring losses, 
but clearly the motoring losses are not the same as the losses under firing 
conditions. Some of the reasons are the lower pressure acting on piston rings 
and bearings, the lower temperatures of the piston and cylinder bore surfaces 
and thus the greater oil viscosity, the greater running clearance of the 
piston, and the missing exhaust blowdown period. 

Nevertheless, a breakdown analysis of motoring losses supplemented by 
experiments on piston and ring friction rigs can be used to identify the 
relative importance of the many components of the total friction and their 
response to changes in design variables. In general, the components of the 
losses expressed in terms of mep tend to fall into three groups: 

(i) Losses due to boundary lubrication, approximately invariable with 
speed. There losses are undoubtedly influenced by the compression 
ratio. 

(ii) Losses associated with oil film friction, varying directly as speed. 
All major rotating parts fall into this group. 

(iii) Losses associated with air pumping, varying as the square of the 
speed. 
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Therefore, the motoring losses can be expressed in the form 
F = A + BN + CN 2 (5-39) 

where F are the losses in mep, N is engine speed, and A, B and C are 
constants. 

Millington and Hartles [20] have measured motoring losses on a large 
variety of automotive diesel engines during the course of development of 
prototype engines. Their work suggests a readjustment of equation (5—39 ) 
coupled with suitable selection of the constants as follows: 


F = A + 


7 -° 1000 


1 ’ 5 ( 1000 ) 


(5-40) 


where F = motoring losses, psi mep 

A = compression ratio minus 4 for a DI diesel 
N = engine speed, rpm 
V = mean piston speed, fpra 

Equation (5~40) represents a sound empirical correlation of the motoring 
loss data obtained from diesel engines. Therefore, it will be used to obtain 
brake quantities from the indicated quantities computed in the engine 
simulation. 
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6. BASIC ASSUMPTIONS OF OTHER COMPONENT MODELS 



7. The wastegate bleeds off a fraction of the exhaust gas stream (rh ), such 

w 

that the flow bypasses the turbocharger turbine but not the power turbine. 

8. The wastegate valve behaves as a throttle valve of variable effective 
throat area. The size of the opening can be governed by the intake 
manifold pressure, p . 
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9. The EGR flow is either specified as a fraction of the intake flow, or the 
EGR valve is modelled as a throttle valve of fixed effective throat area. 

10. The flows through the wastegate and EGR valves are subsonic and adiabatic, 
and the throat pressures are equal to the downstream pressures. 

11. The two exhaust gas streams, through and by-passing the turbine, (m and 

I* 

m ) reunite with no loss of thermal energy or pressure. 

12. The inside wall temperature of all components in the system are known. 

13- There are pressure drops between the following: 

compressor discharge and intercooler inlet (p^ & p^) 

intercooler inlet and intake manifold (p_, & p_) 

3 b 

exhaust manifold and turbine inlet (p^ & pg) 

turbine discharge and power turbine inlet (p^ & p^) 

power turbine discharge and atmosphere (p & p ) 

1 * asm 

14. Wastegate upstream pressure equals turbine inlet pressure. 

15. Wastegate downstream pressure equals turbine discharge pressure. 

16. Turbine and compressor speeds are equal. 

17- The turbocharger (turbine, compressor and shaft) has fixed values of 
rotational inertia and damping, relative to its housing (J and B, 
respectively). 

18. The power turbine shaft is connected to the engine crankshaft via a 

specified gear ratio transmission, thus the power turbine speed is always 
determined directly from the engine speed. 
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7. MODELLING OF SYSTEM COMPONENTS 


7.1 State Variables 

The turbocharged diesel engine system described in Sec. 6 can be modelled 
mathematically by a set of simultaneous, non-linear differential equations 
involving a minimum set of dependent variables. These variables, called the 
state variables, are sufficient to describe completely the system at any point 
in time. Although their number is fixed by the system, there are many choices 
of state variables, some mere convenient than others. The mass, temperature, 
absolute pressure and fuel fraction states of each manifold and the 
turbocharger speed were chosen as the nine state variables of this system. 
Crank-angle was chosen as the independent variable, rather than time. Since 
we are studying constant speed performance, time and crank-angle are directly 
proportional. All mass flow rates (m) in the following analysis refer to mass 
flow per unit crank angle. 

7.2 Derivation of Differential Equations 
7.2.1 Manifolds 

The four state equations for each manifold are derived from four 
interdependent relations for the manifold control volume: 

1. conservation of total mass 

2. conservation of fuel mass 

3. conservation of energy (first law of thermodynamics) 

4. ideal gas law 

The general manifold control volume model is shown in Fig. 4. The total 
mass and fuel fraction state equations for the manifold come directly from the 
continuity relations derived earlier (Equations (3-1) and (3-5), 
respectively) : 
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( 7 - 1 ) 


,dni' 

V 

m 


E m. 

J J 


(~) 

at 


m . 

E - J (F.- F ) 
. m 1 m 
J m 


(7-2) 


where the subscript m refers to the manifold control volume, and the subscript 
j refers to the j-th mass flow that enters or leaves the control volume. 

For the intake manifold, the mass flows entering and leaving the manifold 
include: 


i) compressor mass flow (m c in Fig.1) (entering) which is found from the 
compressor map 

ii) engine intake mass flow (m.^ in Fig. 1) (leaving) which is found 
from the engine model 

iii) EGR mass flow (ra^ in Fig. 1) (entering) which is found frcm the 
EGR valve equations 

For the exhaust manifold, the mass flows entering and leaving the 
manifold include: 

i) turbine mass flow (m^. in Fig.1) (leaving) which is found frcm the 
turbine map 

ii) engine exhaust mass flow (m in Fig. 1) (entering) which is found 
from the engine model 

iii) EGR mass flow (- m^^ in Fig. 1) (leaving) which has the opposite 
sign to intake EGR flow) 
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iv) wastegate mass flow (m in Fig. 1) (leaving) which is found from the 

w 


wastegate valve model 

The state equation for temperature is derived from the general temperature 
state equation (3-26) applied to the manifold control volume. 


(--) , 1 [ft 

v dt' A LV dt' 
m m 


(B - 


h ) 
m 


* 

m 


1_ 

m 


( l m h - Q)] 
m j J J 


(7-3) 


where A, B, and C are defined in (3~27), (3~24), (3~28) and are calculated by 
the thermodynamic property subroutine, 


Q = heat transfer rate to the manifold walls 
and d4>/dt is related to dF/dt and F by ( 3 - 8 ). 

The state equation for pressure is derived from the general pressure 
state equation applied to the manifold control volume (Eq. 3~21 ) : 


_ — £ r_ 1 ^£ (^1) + 1 _ 1 ^£ (-i)i (7-4) 

W 3p/3p L p 3T ^dt' m ^dt ; p 3<j> dt' J w ' 

m m m r m 

The mass, fuel fraction and temperature derivatives are calculated in 

equations (7-1), ( 7" 2) , and ( 7— 3 ) ; the equivalence ratio derivative, d<(>/dt, is 

related to the fuel fraction derivative by Eq. (3-8) ; and the density 

derivatives are calculated by the thermodynamic property subroutine. 

7.2.2 Turbocharger 

The state equation describing the turbocharger is derived from the 
conservation of energy relation. 

--(E . ) = W + W. . . (7-5) 

dt t/c compressor turbine 

where "(E . ) = Jw ~ + Bw 2 (7-6) 

dt t/c dt 

E ^ = total mechanical energy of turbocharger rotor 
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J = rotational inertia of turbocharger 

B = rotational damping of turbocharger 

oi = angular velocity 


W = m (h - h ) 

compressor c 1 2 


Vi , . = nr (h Q - h n ) 

turb me t 8 9 


Solving for doi/dt gives: 


(7-7) 

(7-8) 


dt ' k l( ”o (h 1 ~ V * \ (h 8 - V - B “ 2 ' W-9> 
The enthalpy changes across the compressor and turbine are calculated using 
the compressible flow relations. 


. R/c 

h, P 2 P 

<VV ■ TV 

c H 


R/c 

p 8 P 

(hg hg) - n fc hg [( p ^) 1] 


( 7 - 10 ) 


(7-11) 


The complete state equation for turbocharger is found by combining equations 


(7-9), (7-10), and (7-11). 


7.3 Wastegate and EGR Valve Models 

The wastegate and exhaust gas recirculation (EGR) valves are modelled as 
adiabatic throttle valves of known effective throat area (which is the product 
of actual cross-sectional area and a discharge coefficient. The mass flow 
rates through the valves are calculated using (5-1) and (5-2). 

To determine the EGR flow rate, either the percentage of the engine inlet 
flow which is EGR is specified, or the EGR valve opening is treated as a fixed 
effective-area orifice. 
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The wastegate valve opening area varies between zero and a fixed maximum 
value, and is determined by an actuator linked to a pressure sensor. 

The sensor measures the gage pressure of the intake manifold. At present, a 
simple actuation model has been adopted which gives a continuous closed-form 
relation between effective valve throat area and intake manifold pressure. 

This model can be replaced when the actual actuator characteristics are 
available. 


7.4 Intercooler Model 

The intercooler, which is situated between the compressor discharge and 
intake manifold, serves to increase the density of the charge air by lowering 
its temperature. The intercooler is modelled as a heat exchanger of fixed 
area, heat-tranfer coefficient and cooling flow rate. The change in charge 
air temperature is determined from the non-dimensional heat exchanger 
effectiveness, e. 


e = 


(T - T ) 
__h1_ _h2_ 

(T ~- T ) 

^ hi cr 


(7-12) 


where h refers to the charge air flow (hot) 
c refers to the coolant flow (cold) 

1 and 2 refer to inlet and exit conditions respectively 
T h1 = compressor discharge temperature, T^ 

T^ = intercooler discharge temperature, T^ 

T ^ = coolant inlet temperature (assumed to be fixed) 

Heat exchanger effectiveness is either known as a design parameter or can 
be derived from graphical correlations that are available for the various 
typical heat exchanger configurations. From the latter e can be determined as 
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a function of capacity rate ratio, C . /C , and number of heat transfer 

min max 

units, N , where 

C min = the lesser of the two values of (assumed to belong to the 
charge air flow), 

C = the greater of the two values of me (assumed to belong to the 

ITfciX p 

coolant flow and to be fixed) , 


N = AU/C . 

A min 


A = heat exchange surface area (fixed) 


U = heat transfer coefficient based on A 


Fig. 5 shows a graphical correlation for a cross-flow heat-exchanger with 

unmixed flows, with various ratios of C . /C . Currently the assumption is 

min max r 

being made that C is much larger than C . . With this assumption the 
tifax min 

expression for effectiveness reduces to the following simple form: 


e = 1 - expC-N^) 


(7-13) 


7.5 Manifold Heat Transfer 

With the incorporation of an intercooler, heat transfer frcm the intake 

manifold to the environment becomes small enough to be neglected. Heat 

/ 

transfer from the exhaust manifold, however, cannot be neglected. Heat 
transfer from the gas in the exhaust manifold to the environment involves a 
combination of forced convective 1 heat transfer from the gas to the inside 
walls, conduction from the inside walls to the outside walls and to the water 
jacket, and natural (and probably forced) convection from the outer surfaces 
to the environment. 

It is clearly a complicated process to determine fully the heat transfer 
from inside gas to outside air, involving many detailed assumptions about 
engine geometry and environmental flow conditions. To simplify this task in 
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the current generation simulation program we will deal only with the modelling 
of the heat transfer from manifold gas to interior walls, which necessitates 
that a value for the wall temperature be assumed. 

The heat transfer coefficient is derived from an experimental correlation 
that relates Nusselt, Reynolds and Prandtl numbers for fully-developed 
turbulent flow in circular tubes with large temperature gradients (Eq. 7-14). 


UD 


0.8 0.4 

VDp c u 
L ' ' £L 'i 


Nu . ^ = 0.035 (--) (-*-) 

b f b 


where D = exhaust manifold diameter 

, V = average velocity of gas 

The gas properties are: 

k = thermal conductivity 
u = absolute viscosity 
p = density 


0 ^ = specific heat at constant pressure 
Subscripts are: 


b: refers to properties at T . 
f: refers to properties at T f . 
Temperatures are: 

T fa = bulk (averse) temperature 
T w = wall temperature 
T^ = "film" temperature, 


T = -(T 
f 2 w w 


+ T b } 


(7-14) 


(7-15) 
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7.6 Pressure Losses 

Pressure loss terms have been included at five locations in the overall 
system model. These are: 

between compressor discharge and intercooler inlet, 
across intercooler, 

between exhaust manifold and turbine inlet, 
between turbine outlet and power turbine inlet 
between power turbine outlet and atmosphere. 

Each of these pressure drops is calculated using the corresponding friction 
factors and friction coefficients for the geometry of each passage. 

For straight-sections: 

P = 4f(L/D)(pV 2 /2) (7-16) 

where L = length of passage 
D = diameter of passage 
p = bulk density 
V = average velocity 

f = friction factor, which is correlated by (7-1 7) for the surface 
roughness and range of Reynolds numbers to be encountered. 

Note . that 

f = 0.046(Re)^' 2 (7-17) 

with Re = pVD/p 

For bends, enlargements, -contractions, etc: 

P = K fP V 2 /2 (7-18) 

where K f = friction coefficient for a particular passage geometry. Values 
of K^, for typical geometries are commonly available. 
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8. THERMODYNAMIC PROPERTIES 


Our thermodynamic model assumes that the various control volumes contain 
mixtures of air aid combustion products throughout the total engine system. 

By utilizing the concept of the instantaneous average equivalence ratio 
defined in Section 3.1, the contents of any control volume can be represented 
as one continuous medium. Furthermore, assuming ideal gas behavior and 
thermodynamic equilibrium, the instantaneous gas properties can be determined 
from a knowledge of pressure, temperature and average equivalence ratio in the 
control volume. 

When the temperature of the cylinder contents is below 1000 K, they are 
treated as a homogeneous mixture of non-reacting ideal gases, their properties 
being calculated using the procedure outlined below [21]: 

The hydrocarbon-air combustion reaction is written as: 

£<}> CH +0 + <Jj N -*• x. CO + X-H-0 + x_C0 + x,.H„ + X..CL + t|>N (8-1) 
y I 223 4202 

where i|> = the molar N:0 ratio of the products, 

y = the molar H:C ratio of the fuel, 

<() = the average equivalence ratio, 

x. = moles of species i per mole of 0 o reactant 
1 d 

and e = 4/(4+y) (8-2) 

The quantities x^ are determined by using the following assumptions: 

a) for lean mixtures (<j> £ 1) can be neglected. 

b) for rich mixtures (<p > 1) 0^ can be neglected. 

c) for rich mixtures, the gas water reaction 

co 2 + h 2 -*• co 2 + h 2 o 

is in equilibrium with equilibrium constant K(T). 
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The solution for the x.^ is shown in Table II, where C is obtained by- 
solving equation (8-3) for its positive root. 

(1 - K)C 2 + 2[1 - e<p + K(<f> - 1 + e<j>)]C - 2Ke<p( <p - 1) = 0 (8-3) 

The value of K(T) is obtained by curve fitting JANAF table data over the 
temperature range 400 to 3200 K and is given by 

ln(K(T) ) = 2.743 - 1 .761 /t - 1.611/t 2 + .2803/t 3 (8-4) 

where t = T/1000, and T is the temperature in Kelvins. 

If the grams of products per mole of 0 reactant is expressed as 
M = (8e + 4)4* + 32 + 284; (8-5) 

the specific enthalpy h and the specific heats at constant pressure and 


constant composition, c^ and c^ respectively, can be expressed by the 


following relationships: 
6 4 


1 ' pJ a p5 

h = - Z x. £ (a. . t - -~ + a.,) 

M i= i 1 j=1 ij J t 16 


( 8 - 6 ) 


1 b 4 . ■ . a. c 

c = - Z x. z (a. ,t J_ + -~) 
p M >1 1 >1 « t 2 


(8-7) 


. 6 9x. 4 j a 

c, = - 1 r-- Z (a. . t - + a.,) - 

<t> M . =1 94» ij j t 16 


1 *• E (a. . ~~ + a. fi ) 

M 2 i=1 1 j.i 1J J 16 


( 8 - 8 ) 


The coefficients a„ are obtained by curve fitting JANAF table data to the 
above functional form. The values of a_ are given in Table III. The 


resultant c is in cal/g-K, while h and c, are in kcal/g. 
P 4> 


Since the cylinder contents are being treated as a mixture of non- 
reacting ideal gases, the density of the mixture is given by 
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(8-9) 



where 

R = the universal gas constant (1.9869 cal/mole-K) 

o 

and M, the average molecular weight of the mixture, is given by 
M/((1-e)<f> + 1+iJj) <p £ 1 

and (8-10) 

M/((2 - e)«j> + 4<) <fi > 1 

Then, the partial derivatives of the density with respect to 
temperature, pressure, and equivalence ratio are given by 



( 8 - 11 ) 


Ifi = £ 

3p p 


( 8 - 12 ) 


^£ = i£ f 8-1^1 

d<p ' m d * 

'When the temperature of the cylinder contents is above 1000 K, their 
properties are calculated with allowance for chemical dissociation, according 
to the calculation method described in [22]. This is an approximate method 
based on curve fitting data obtained from detailed thermochemical calculations 
[23] to a functional form obtained from concentrations within the burned gases 
are not calculated, the bulk thermodynamic properties needed for cycle 
analysis are accurately determined. 
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9. TRANSPORT PROPERTIES 


The heat transfer correlations relate the heat transfer coefficient to 
the Reynolds and Prandtl numbers and the thermal conductivity. The 
calculation of the heat transfer rates will therefore require values for the 
viscosity and the Prandtl number (from which the thermal conductivity can be 
obtained). We have used the approximate correlations for the viscosity and 
the Prandtl number of hydrocarbon-air combustion products developed by 
Mansouri and Heywood [24]. 

The NASA equilibrium program [23] was used to compute the viscosity of 
hydrocarbon-air combustion products as a function of temperature, T, 
equivalence ratio, <j>, and pressure p. It was shown that the viscosity of the 
combustion products was satisfactorily correlated by a power-law based on air 
viscosity data, corrected for the effect of equivalence ratio, i.e., 

U prod [kg/ms] = 3.3x1 0 _7 T°* 7 /(1 + 0.027<j>) 

(9-1) 

for 500 K i T S 4000 K 
0 < <j> $ 4 

Note that the viscosity of the combustion products is independent of the 
pressure. 

The equilibrium Prandtl number of hydrocarbon-air combustion products was 
also calculated over the above ranges of temperature, pressure, and 
equivalence ratio. Using a second order polynomial of Y to curve fit the 
above data, it was shown that the following correlation for lean ($ < 1) 
mixtures predicted values in good agreement (within 5%) with the data, i.e., 


45 



(9-2a) 


Pr = 0.05 + 4.2(Y - 1) - 6.7(Y - I) 2 
for 500 K < T £ 4000 K 
and <j> £ 1 

For rich mixtures (<j» > 1), a reasonable fit (less than 10$ error) to the 
equilibrium Prandtl number values calculated with the NASA program was found 
to be the following: 

Pr = [0.05 + 4.2CY - 1) - 6.7(Y - 1) 2 ]/[1 + 0.01 5x10~ 6 (<J>T) 2 ] 

for 2000 K < T £ 3500 K (9~2b) 

and 1 < $ £ 4 
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10. ITERATION PROCEDURE 


1 0.1 Basic Method of Solution 

When the individual submodels of the total engine system are brought 
together to form a complete model, the result is a set of simultaneous 
first-order ordinary differential-equations. To perform predictive 
calculations with the cycle simulation, these equations must be integrated 
simultaneously over the full operating cycle. Note, however, that some of the 
governing equations like the mass flow rate through the intake or the exhaust 
valve and the non-dimensional fuel burning rate apply only during parts of the 
cycle. 

Integration of the governing equations is performed numerically using a 
standardized code developed by Shampine and Gordon [27]. The optimal 
integration time step is calculated automatically in the routine so as to 
maximize efficiency. Detailed documentation of the integration routine is 
provided in the listing of the code. 

10.2 Crank Angle By Crank Angle Integration 

The engine model calculates the state variables in one master cylinder of 
a multi-cylinder engine, while the manifold and other component models have 
inherent multi-cylinder capability. The interaction between the master 
cylinder model and the other components takes place in the manifolds. To 
simulate the effect of additional cylinders on the manifold conditions and 
hence on the entire system behavior, the conditions in the other cylinders are 
assumed to vary as echoes of the first cylinder, shifted by the appropriate 
phase lags. The intake and exhaust mass flow profiles calculated by the 
engine model for the master cylinder are used to represent the profiles for 
the other cylinders. These profiles are updated after the completion of each 
engine cycle for use in the next iteration. 
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A flow chart showing the overall structure of the engine cycle 
simulation is shown in Figure 6. After the initialization of the state 
variables, the simulation proceeds with the simultaneous integration of the 
state equations for the manifolds, the engine and the turbocharger. The 
main program determines which sets of equations the integrator will work 
with, and calculates and prints our results as the cycle proceeds. Four 
secondary routines corresponding to the intake, compressor, combustion, and 
exhaust processes are called in turn by the integration routine to supply it 
with the derivatives of the variables which it must integrate. A third 
level of utility routines, such as the other component state equations, 
thermodynamic and transport property routines, tables of valve flow areas 
and discharge coefficients, is called by the second level routines to help 
in the evaluation of the necessary derivatives at each step. 

Approximate estimates of all state variables at convergence are 
assumed. In general, more than one iteration is required to model an engine 
operation under steady conditions. The integration continues until the 
system reaches a quasi-steady state, which is defined as the state in which 
the value of each state variable at a particular crank-angle is the same (to 
within a specified interval) as it was at that crank-angle in the previous 
engine cycle. Thus, at quasi-steady state a set of repeating patterns of 
variation, one for each state in the system and each with a period of one 
engine revolution, is achieved. 
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11 . PROGRAM INPUTS AND OUTPUTS 

11.1 Inputs 

To operate the simulation, a set of input parameters must be specified. 
These include fixed system dimensions, initial values of state variables and 
computing tolerances. They are listed and briefly described below. 

11.1.1 Engine operating conditions 

i. fuel type 

ii. mass of fuel injected per cycle 

iii. injection timing 

iv. engine speed (rpm) 

v. constants in the ignition delay correlation 

vi. nominal burning duration 

vii. constants related to the burning rate distributions 

11.1.2 System dimensions and design parameters 

a. Engine parameters 

i. number of cylinders 

ii. cylinder bore and stroke 

iii. clearance voIuhb 

iv. connecting rod length 

v. valve timings (crank angles at which the intake and exhaust 
valves open and close) 

vi. tabulated values for the valve open areas and discharge 
coefficients 

b. Other component dimensions 

i. manifold dimensions 

ii. connecting passage dimensions 
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iii. EGR and wastegate valve dimensions 


iv. wastegate actuator characteristics 

c. Intercooler characteristics 

i. coolant flow heat capacity 

ii. heat exchange surface area 

iii. heat exchanger heat transfer coefficient 

d. Turbomachinery parameters 

i. turbocharger rotational inertia 

ii. turbocharger rotational damping 

iii. power turbine gear ratio 

iv. power turbine transmission efficiency 

v. compressor, turbine and power turbine maps 

(see Sec. 11.2 - Turbocharger and power turbine maps) 

11.1.3 Heat transfer and turbulent flow parameters 

a. constants for the Nusselt-Reynolds number correlation 

b. the temperatures of the piston, cylinder head, cylinder walls and 
manifold walls 

c. the turbulent dissipation constant 

11.1.4 Engine friction parameters 

Empirical constants associated with the calculation of friction losses 

11.1.5 Initial conditions 

a. intake manifold pressure and temperature 

b. exhaust manifold pressure and temperature 

c. turbocharger speed 
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11.1.6 Ambient conditions 

a. intake temperature 

b. intake pressure 

c. final exhaust pressure 

11.1.7 Computational parameters 

a. convergence margins for each state variable 

b. error tolerances for integration of the differential equations 

c. other parameters used in the integration algorithm 

(For a detailed description of these parameters, see the engine 
simulation code.) 

11.2 Turbocharger and Power Turbine Maps 

Maps give the interrelationships among mass flow rate, efficiency, 
pressure ratio and rotor speed for each of the three turbonHchinery 
components: turbocharger compressor and turbine, and power turbine. The maps 
are put into the simulation in tabular form and must therefore be reduced to 
tabular form if initially in graphical form. The tables are interpolated 
during the simulation to find the necessary information. The map variables of 
mass flow rate and rotor speed are corrected by factors relating actual inlet 
conditions to standard conditions. The speed correction factor involves 
inlet temperature, and the mass flow rate correction factor involves inlet 
temperature and pressure. 

The compressor map variables are 1) m , 2) p^/p^ 3) n c and 4) w. 

The mass flow rate is normalized according to Pg and Tg. The rotor speed is 
normalized according to Tg. 

The turbine map variables are 1) m. , 2) p Q /p n , 3) a and 4) a>. 

t o y t 
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The mass flow rate is normalized according to Pg and Tg. The rotor speed is 
normalized according to Tg. 

The power turbine map variables are 1) m , 2) p /p , 3) n , 4) u . 

pt I U II P O 

Mass flow rate is normalized according to p 1Q and T 1Q . The shaft speed, which 

is directly determined from the engine speed and the specified gear ratio, is 
normalized according to T^. 

11.3 Outputs 

Four types of outputs are generated by the cycle simulation: 

11.3.1 Input echo 

A listing of all the input parameters, including some quantities derived 
directly from the given inputs (e.g. engine displacement and compression 
ratio) 

11.3*2 Major crank-angle by crank angle results 

At specified crank-angle intervals, the values of the following state 
variables are returned: 

a. cylinder pressure, temperature, and average equivalence ratio 

b. intake manifold pressure, temperature and average equivalence ratio 

c. exhaust manifold pressure, temperature and average equivalence ratio 

d. turbocharger speed 

In addition, the following other quantities are reported at the same 
intervals: 

a. engine total heat transfer rate 

b. heat transfer rates to different surfaces of interest 

c. engine work done 

d. turbulent flow model results (such as mean flow velocity, turbulent 
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d. turbulent flow model results (such as mean flow velocity, turbulent 
Intensity, macroscale of turbulence) 

e. a code which monitors the performance of the integration routine 
Integrating through the different processes for the master cylinder, the 

following quantities are reported during the corresponding process: 

a. Intake 

i. mass flows through each valve 

ii. mass fraction of fresh charge (air) in the cylinder 

b. Combustion 

i. a non-dimensional fuel burning rate 

ii. fuel burnt as a function of total fuel injected 

c. Exhaust 

i. mass flow rate and velocity through exhaust valve 
11.3-3 Integrated results and cycle performance 

After completion of an engine cycle, a summary of results obtained by 
integrating some of the governing equations over the cycle is given. 
Integrated results include the following: 

a. volumetric and thermal efficiencies 

b. gross indicated, pumping and friction mean effective pressures 

c. total heat loss 

d. ignition delay period 

e. estimated mean exhaust temperature 

f. mass in cylinder at IVO and at IVC 

g. mass of air inducted per cycle 

h. total heat and work transferred during each process 

i. results of an overall energy balance 
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11.3.4 Sub-npdel results 

After the overall cycle results are listed, detailed results for the main 
sub-models of the cycle simulation are given at specified crank angle 
intervals. These quantities include the following: 

a. total engine intake and exhaust mass flow rates 

b. compressor, turbocharger turbine and power turbine mass flow rates, 
pressure ratios and efficiencies 

c. wastegate (turbocharger turbine bypass) and EGR mass flow rates 

d. power turbine work transfer 

e. pressures and temperatures at various system locations 

f. intake and exhaust manifold heat transfer rates 

g. intercooler effectiveness 
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Table I: Constants for Arrhenius Equation for Ignition Delay [ms] 



Cetane 



E/R 



Fuel 

No. 

A 

n 

kelvins 

Ref. 

Notes 

diesel 


53.5 

1 .23 

676.5 

(a) 


n-cetane 

100 

0.872 

1.24 

4050 

(a) 

Bomb data at 
high pressure 

n-heptane 


0.748 

1.44 

5270 

(a) 


diesel 

45-50 

4.05x1 0 -2 

0.757 

5470 

(b) 

Steady flow 

diesel 


2.43x1 0' 6 

2 

20915.3 

(c) 

reactor 

kerosene 


1 .68x1 0 ^ 

2 

19008.7 

(c) 

Continuous 

cetane 


4.04x1 0 -10 

2 

25383 

(c) 

flow 

apparatus 

diesel 


3.45 

1 .02 

2100 

(d) 


(a) Igura, 

S. , Kadota, T. , Hirozasu, H., 

"Spontaneous 

Ignition Delay of Fuel 


Sprays in High Pressure Gaseous Environment," Jap. Soc. Mech. Engrs., Vol. 
41, No. 345, p. 1559, (1975). 

(b) Stringer, F.W. , Clarke, A.E., Clarke, J.S., "The Spontaneous Ignition of 
Hydrocarbon Fuels in a Flowing System," Proe. Auto. Div. , Institution of 
Mech. Eng., 1970. 

(c) Spadaccini, L.J., and TeVelde, J.A., "Autoignition Characteristics of 
Aircraft-Type Fuels," United Technologies Research Center, NASA report 
CR-1 59886, 1980. 

(d) Watson, N. , Pilley, A.D. and Marzouk, M. , "A Combustion Correlation for 
Diesel Engine Simulation," SAE paper 800029, 1980. 
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Table II: 


Burned Gas Composition under 1000°K 


x. (moles/noleO_ reactant) 
1 2 . 


■ 

Species 

4> £ 1 

$ > 1 

■ 

co 9 


£$ -C 

2 

H,0 

2 (!-£)$ 

2 (l-e£) + c 

3 

CO 

0 

C 

B 

H 2 

0 

2W-1) - C 

5 

°2 

1 - $ 

0 

6 

N 2 

* 

'p 


sum 

(l-£)d + 1 + $ 

(2-£)<t> + P 
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Tabic III: Coefficients for Polynomial Fit to Thermo dynamic Properties 


Coefficients for ICO 0 -: < T < 500°K 


■ 

Species 

3 il 

a i 2 

S i3 

a ii 

a . _ 
10 

a. * 
16 

1 

CO, 

IB 

1 

-11.232 

m 

. 006 7 6 7 

-93.75793 

9 

H ,0 

7 . 509 7 

B 

3.41S7 

-1.1790 

.001436 

-57.08004 

3 

co - 

6.9739 

- . S23S3 

2.9420 


.0004132 

-27.19597 

n 

H 2 


.16170 

- . 21S21 

?0 ftO-7 

• s U U m 

-.016252 

-.113189 

5 

°2 

mm 

mm 

j||||||j 

-.32674 

.004359 

.103637 

6 

\T 

7.0922 

-I.295S 

3.2069 

- 1.2022 

-.0003455 

-.013967 


Coefficients for 5C0°K < T < ^OGO°K 


■ 

Species 

a il 

a i 2 

a i3 

a i4 

a i5 

2 ^ 

a i 6 


C0 2 

11.940 

2.03S6 

-.47029 

.037363 

- .58945 

EB 

2 

« o 
H 2 ° 

6.1391 - 

4.6073 


.066695 


->u • l; o o -] 

3 

CO 

7.0996 

B 

- . 2S775 

.022356 

-.15937 

-27.73464 

D 

H, 

H2SI 

Qjg 

-.23313 

.019515 

.16113 

. 7649S 


°2 

7.8658 

.63537 

-.031944 

-.0026370 

-.20139 

-.393455 

6 

N 2 


Ba 

- . 32S99 

.025610 

-.11395 

- . 331S35 


* picked to give enthalpy datum at Q°K 
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Figure 1: Turbo-compounded diesel system configuration 
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Figure 2: Typical rate of heat release diagram 

for a DI diesel engine 
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Figure 3: Energy cascade model 
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Figure 4: General manifold control volume 
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Figure 6: Flowchart of simulation program 
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